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We report constant- volume and constant-pressure simulations of the thermodynamic and dynamic 
properties of the low-temperature liquid and crystalline phases of the modified Stillinger- Weber 
(mSW) model. We have found an approximately linear increase of the effective Gaussian width of 
the distribution of inherent structures. This effect comes from non-Gaussianity of the landscape and 
is consistent with the predictions of the Gaussian excitations model representing the thermodynamics 
of the configurational manifold as an ensemble of excitations, each carrying an excitation entropy. 
The mSW model provides us with both the configurational and excess entropies, with the difference 
mostly attributed to vibrational anharmonicity. We therefore can address the distinction between 
the excess thermodynamic quantities often used in the Adam-Gibbs (AG) equation. We find a 
new break in the slope of the constant pressure AG plot when the excess entropy is used in the AG 
equation. The simulation diffusivity data are equally well fitted by applying a new equation, derived 
within the Gaussian excitations model, that emphasizes enthalpy over entropy as the thermodynamic 
control variable for transport in viscous liquids. 
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I. INTRODUCTION 

Both thermodynamic and dynamic properties of 
structural glass-formers are unusual and not fully 
understoodiii^ It has long been suggested that puzzles 
of the dynamics of supercooled liquids can be unrav- 
eled from the properties of their energy landscapes;^, At- 
tempts to build a description of dynamics based on the 
liquid thermodynamics go back to the Adam-Gibbs (AG) 
theory^ which emphasizes configurational entropy as the 
origin of non-Arrhenius dynamics. Even though the en- 
tropic paradigm seems to be currently most successful 
in describing relaxation^^"?*'^ other approaches have em- 
phasized the energetic aspects of the problem in terms of 
the activated kinetics over enthalpic barriers increasing 
in height with lowering temperaturei^i^ii^iiSiiii Some very 
recent datai^ seem to disfavour the existence of diverging 
lengthscale assumed in "cooperative region" models-'^ 
supporting instead the picture of activation-based dy- 
namics in glass-forming materials. 

Configurational entropy is a property not directly 
accessible from laboratory experiment, but is increas- 
ingly available from computer simulations of model 



systems 
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Excess entropy, s^y(r), i.e. the 



liquid entropy over that of the thermodynamically stable 
crystal, has been successfully used^^'^" instead of config- 
urational entropy Sp y (T) in the AG relation 



ln(Tp,y(T)/ro) - 



Ts%y{Ty 



(1) 



where Tpy(T) is the time of structural relaxation and 
To is the time characteristic of quasi-lattice vibrations. 
The subscripts P, V in the entropy and relaxation time 
specify the conditions, constant-pressure (P) or constant- 



volume (V), at which the entropy and the relaxation time 
are measured. Historically, the definition of the config- 
urational entropy used by Adam and Gibbs^ was that 
of the full configurational entropy including the entropy 
of basin vibrations. However, since they explicitly de- 
manded the vibrational entropy to cancel between the 
liquid and the crystal, their configurational entropy is in 
fact the entropy of inherent structure introduced by Still- 
inger and Weber (see below), We will therefore refer to 
the configurational entropy in this latter definition. 

The full configurational entropy and Stillinger's en- 
tropy of inherent structures are in fact not equiv- 
alent. Most studies, both laboratory22iS^i2i and 
computational,"i2ii£dS, have shown that the excess en- 
tropy has a significant contribution from the vibra- 
tional manifold related to the excess density of vibra- 
tional states in liquid, g^Py{oj), compared to the crystal, 

gp^^{uj). The excess vibrational entropy from harmonic 
motions is a sum over the vibrational spectrum: 



So long as the structure does not change, the ex- 
cess density of states is independent of temperature for 
purely harmonic vibrations resulting in a temperature- 
independent excess harmonic entropy Aspy and zero 
contribution to the excess heat capacity. The anhar- 
monicity of atomic and molecular vibrations leads to two 
effects: (i) sample expansion at constant-pressure heat- 
ing and (ii) deviations of the vibrational excess entropy 
Asp y from the harmonic formula in Eq. The first ef- 
fect makes the density of quasi-harmonic vibrations vary 
with temperature and the second effect requires extract- 
ing the thermodynamics of vibrations without relying on 
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FIG. 1: Excess (for laboratory data) and configurational 
(for computer simulations) entropy (a) and heat capac- 
ity (b) of some glass-forming liquids: glass-forming alloy 
Pd43Cu27NiioP20 (Pd43)r^ o-terphenyl^S, (OTP) and their 
computer models, binary mixture Lennard- Jones (BMLJ)^ 
and Lewis and Wahnstrom o-terphenyl (LW-OTP).— The 
closed points and thick lines are laboratory data, the open 
points and dashed lines are computer simulations, and the 
thin lines are fits to the IG excitations model (from Ref. H^)- 
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FIG. 2; Adam-Gibbs plot of the experimental dielectric re- 
laxation data of oterphenyl (OTP)^ and simulated self- 
diffusivity^^ of Lewis and Wahnstrom o-terphenyl (LW-OTP, 
NVT simulations). Excess entropy (at constant pressure) was 
used for experimental points and configurational entropy (at 
constant volume) for simulated data. The dashed line is a 
linear fit through the data with the intercept at tq — 10"^** s. 
The simulation data are shifted vertically to cross the vertical 
axis at about the same relaxation time. The Adam-Gibbs re- 
lation does not linearize the dielectric data in the whole range 
of temperatures displayed,^" and the plot only serves to illus- 
trate the difference in the range of Adam-Gibbs parameter 
{Tsp''y{T))~^ and relaxation times explored by simulations 
{V = Const) and experiment (P — Const). 



the harmonic approximation [Eq. 

For most systems studied to date, the density of quasi- 
harmonic vibrations is known to depend weakly on tem- 
perature except for the low-frequency feature known as 
the Boson peak. The excess vibrational density of states 
from this region in liquid selenium was shown to produce 
ca. 30% of both the excess entropy and heat capacityi^ 
These results, and the puzzling ability of both the config- 
urational (simulations) and excess (laboratory) entropies 
to fit the relaxation data according to Eq. ([1]), have 
lead to the suggestion that configurational and excess 
entropies of glass formers might be proportional to each 
other.25'26 There are simulation^^ and laboratory^^ data 
in support of this proposal, although the question is still 
not fully settled^S 

The resolution of the problem of partitioning the ex- 
cess entropy and heat capacity between vibrational and 
configurational manifolds is not straightforward. Re- 
cent experimental data have argued in favor of both a 
negligible contribution of vibrations to the excess heat 
capacitj^ and a significant contribution amounting from 
30% (Refs. P and ^ to 50-60% (Refs. '33 and 3?). 
Computer models do not tend to help much in resolving 
this issue due to a general disconnect between simulated 
and experimental results. The situation is illustrated in 
Fig. [1] where excess (laboratory experiment) and config- 
urational (computer simulations) thermodynamic quan- 
tities are compared for two substances, o-terphenyl and 
the glass-forming alloy Pd43Cu27NiioP2o, for which both 
force-field models^i^^ and experimental result a^^i'^° are 



available. For o-terphenyl, there is a major difference in 
both the magnitude and temperature dependence of the 
entropy and heat capacity, though this is not unexpected 
given that a flexible 14-carbon molecule is being mod- 
eled with a rigid three-bead particle (making it much 
less entropic). The difference is not removed by com- 
paring the experiments with the limited simulation data 
available for constant pressurei^ For Pd43Cu27NiioP2o 
glass-former, as modeled by Kob and co-workers binary 
mixturCfSi the agreement is better and a noticeable dif- 
ference exists only for the heat capacities, however the 
need for many components in the experimental analog is 
a disadvantage. 

Clearly, to resolve the question of vibrational contri- 
butions to the excess heat capacity, the field is in need of 
better model systems. A minimum need is for a system in 
which both excess and configurational data are available, 
something so far lacking in computer models (except for 
water— which is famously anomalous and unrepresenta- 
tive of the glass problem we are addressing) . The recent 
exploration of the Stillinger- Weber type model by Mo- 
linero et al^ offers an opportunity to fill this gap since 
the crystal state is available and, for chosen parametriza- 
tions, the liquid can be supercooled without limit. Our 
simulations here take advantage of this opportunity. 

Many recent simulations have offered access to both 
configurational thermodynamics and dynamics reporting 
the success of the AG relationi^ii^ii^i^l [Eq. As men- 
tioned above, the low-temperature part of experimental 
relaxation data can be fitted by the same relation where 
s'^py{T) is used instead of Spy(T).^° While this creates 
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a puzzling contradiction, one needs to keep in mind the 
difference in the scales of the two sets of data, which 
is illustrated in Fig. [5] comparing the laboratory dielec- 
tric datai^S for OTP with the results of simulationslS' for 
the Lewis and Wahnstrom (LW) model of OTP,^ The 
difference in slopes of simulations and laboratory data 
may originate from the use of different ensembles, con- 
stant volume and constant pressure, respectively, and of 
different entropies, configurational and excess. The AG 
relation in fact linearizes the laboratory data only at low 
temperatures, below the crossover temperature asso- 
ciated with either the mode-coupling critical temperature 
or the temperature Tf, of the Stickel analysis. The ex- 
perimental data can be linearized in the AG plot with 
different slopes below and above Tx, as we also describe 
below for the modified Stillinger- Weber (mSW) model, 
when the excess entropy is used in the AG plot. 

Both the range of temperatures and the property stud- 
ied can affect the conclusions regarding the validity of 
the AG relation from the simulation data. This is il- 
lustrated for SPC/E water in Fig. [3l The dielectric re- 
laxation times collected from Molecular Dynamics (MD) 
simulations in a broad range of temperatures^^ show a 
break in the slope of the AG plot (configurational en- 
tropy is taken from Ref. [27l ). The change in slope is 
much less pronounced for the one-particle rotational and 
translational relaxation times (triangles and diamonds 
in Fig. [3]) than for the many-particle Debye relaxation 
time (circles in Fig. [3]) . Diffusivity was the only dynamic 
property considered in the original report of success of 
the AG relation for SPC/E water pi which indeed shows 
almost linear AG plot. Notice that the change of the 
slope seen in the simulation data in Fig. [3] is uniformly 
observed in laboratory dielectric measurements of a num- 
ber of glass-formers.^- However, a downward break in the 
slope, in contrast to the upward-curved dependence seen 
for SPC/E water in Fig. [31 is typically obtained for molec- 
ular glass-formers, as is also the case with the AG plot of 
the mSW model discussed below. In case of water, this 
discussion should be supplemented by the question of the 
relevance of the AG plot to the approach to a thermody- 
namic (near-critical) singularity;^ 

Two of us have recently suggested to describe the ther- 
modynamics and dynamics of glass-forming liquids close 
to Tg in terms of configurational excitations with a Gaus- 
sian distribution of excitation energies (Gaussian excita- 
tions model) (^i^ Entropy, calculated from microcanon- 
ical ensemble, was the starting point for developing the 
theory. Each excitation was assumed to carry the excess 
entropy sq which was subdivided into configurational and 
vibrational parts, sq = s^^^ + Sq. These excitation en- 
tropies and the number of configurations of an ideal gas 
of excitations provide the excess and configurational com- 
ponents of the entropy. The excess entropy is connected 
to the configurational entropy by the relation 

SpV(T) = s'^py[x{T), T] + x{T)s-^\ (3) 
The effect of the vibrational excitation entropy Sq^^ is not 
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FIG. 3: Diffusivity (triangles), relaxation time of one-particle 
rotational diffusion (diamonds), and Debye relaxation time 
(circles) of SPC/E wate^^a vs l/jTs^/iT)) with configura- 
tional entropies taken from Ref. 1271 . The dotted lines con- 
nect the points, the solid lines are linear regressions through 
high- and low-temperature portions of the dielectric relax- 
ation data. Tm = 273 K indicates the melting temperature of 
laboratory water at ambient conditions. 

limited to the last term in Eq. ([3]) since the excited-state 
population x{T) is determined by the total excitation 
entropy sq and the excitation energy cq: 

x{T) = (l + exp[so-/3(eo-2x(r)A)])-\ (4) 

where f3 is the inverse temperature. The energies of exci- 
tations are assumed to belong to a Gaussian distribution 
with the width 2AT and the Gaussian width parameter 
A. This energy lowers the energy of excitations in Eq. 
(HI) to the extent determined by the excited state pop- 
ulation and A, e(T) = eo — 2x{T)X, thus producing a 
self-consistent equation for x{T). 

The excess vibrational entropy originating in the ex- 
cess density of states of a liquid relative to its crystal 
provides an extra driving force, in addition to the larger 
number of basins, for the system to reach the top of 
its energy landscapci^ This extra entropic driving force 
will lead to an increased fragility of the liquid if the ex- 
cited state population change occurs in the temperature 
range of interest (near Tg). In such a case, changes in 
the excess density of vibrational states will be observ- 
able directly through neutron scattering study of glasses 
of different Active temperatures, as has been reported 
for some cases. However, there are strong theoretical 
suggestionsiii^^ that in the case of very fragile liquids 
the excited state is highly populated above Tg (due, in 
principle, to an earlier phase transition below T^^i^) in 
which case the excess entropy due to vibrations (which 
can be quite significan t^^'^^ ) will change little with tem- 
perature. There will be a vibrational entropy difference 
from the crystal due to different vibrational density of 
states [Eq. but there will be no temperature depen- 
dence to this excess. The higher heat capacity of the 
fragile liquid must then have some other source. One 
such source could lie in the T-dependence of the Gaussian 
width of the excitation profile, which will be addressed 
below. 
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This conclusion is, of course, limited by the assump- 
tion that Sq'*^ is temperature- independent. This assump- 
tion in fact implies the neglect of vibration anharmonic- 
ity making eigenvalues of the landscape basins depend 
on temperature, also including the effect of thermal ex- 
pansion. If this assumption is violatedfS a vibrational 
contribution appears in the excess heat capacity. We 
found, however, that the assumption of Sg'*^ — Const 
describes laboratory fragile liquids well and also found 
that common model fluids used in simulations classify as 
intermediate/strong in terms of their fragility j^^'^^ 

For intermediate and strong liquids, the tempera- 
ture change in the configurational entropy is driven by 
the changing population of the configurationally excited 
state, as was the case in the original two-state Angell- 
Wong model."'" In this case, the temperature dependence 
of x{T) in Eq. ([3]) contributes to the heat capacity, and 
there is a non-zero vibrational component in the ex- 
cess heat capacity.^^ A nonzero Sg'*^ is therefore sufh- 
cient to produce a vibrational excess heat capacity for 
intermediate/strong liquids, while this property should 
be temperature-dependent for a vibrational excess heat 
capacity to exist for fragile liquid. 

The distinction between the thermodynamics of frag- 
ile and strong glass-formers is attributed in the Gaussian 
excitations model to two parameters: critical excitation 
entropy and critical temperature. In order for the fragile 
behavior to be realized, the excitation entropy should be 
higher than the critical value sqc — 2 and the tempera- 
ture should be lower than Tc ~ A/2. The binary mixture 
Lennard-Jones (BMLJ) liquids, which we have analyzed 
previously;^ and the mSW liquid analyzed here all have 
excitation entropies falling below the critical value and 
thus classify as strong/intermediate liquids. The excess 
heat capacity of the mSW liquid has, therefore, a vibra- 
tional component ~ 15 % as discussed below. 

The dynamic part of the excitation modcP^ has offered 
an expression for the relaxation time alternative to the 
AG formula in terms of the configurational heat capacity 
instead of the configurational entropy: 

DT' 

Hr/ro) = Y-fT^Y 

where D and T' are model parameters varied in fitting 
the experiment. This relation gives the fits of relaxation 
times from laboratory and computer experiments com- 
parable to those based on the AG relation. However, it 
carries a problem similar to the one encountered in ap- 
plications of the AG equation. The experimental data 
can be fitted by using the excess heat capacity, c^?, while 
simulation data can be accounted for with the configu- 
rational part Cp. Here we offer some insights into this 
problem from the data collected for the mSW modeli^ 
Some predictions of the excitations model deviate 
from popular models of landscape thermodynamics. The 
concept of potential energy landscape, proposed by 
Goldstein^ was formalized by Stillinger and Webei™ in 
terms of the thermodynamics of inherent structures re- 



ducing the many-body problem of describing liquids to 
a single "reaction coordinate" </> defined as the depth of 
the potential energy minimum (per molecule, (f> = ^/N, 
N is the number of molecules). The probability to find 
a minimum of depth is determined by the number of 
minima of given depth ^{(f>) and two free energies, the to- 
tal thermodynamic free energy /(T) and the free energy 
f^{4>,T) of the system exploring the phase space within 
the basin surrounding the minimum of depth (/): 

In T)] = -(3cp + u{4>) + /? [/(T) - /"C^, T)] . 

(6) 

Here, w(0) — N ^ ln[ri((/))] is the enumeration function. 
The subscript specifying ensemble (constant volume or 
constant pressure) is omitted in Eq. ^ for brevity. In 
case of constant- volume conditions, f{T) is the common 
notation for the free (Helmholtz) energy per particle. 
At constant pressure, /(T) should be understood as the 
Gibbs energy and (/) is the potential enthalpy minimum4i 
We will drop the ensemble specification in the landscape 
variables below reserving it to the equilibrium properties, 
e.g. Cpy(r), where this distinction is critical. 

The formalism of inherent structures is particularly 
convenient when /''{(/), T) is independent of (/)■ Other- 
wise, Eq. ^ is formally a definition of the basin free 
energy. It was found that at high temperatures acces- 
sible to simulations the harmonic part of /''((/), T) is a 
weak linear function of (j)r^^^ The main focus of the for- 
malism is, however, on the enumeration function. It is a 
decreasing function with lowering 4> allowing the liquid to 
explore higher energy states with increasing temperature 
(entropy drive). It was suggested that Derrida's Gaus- 
sian model^ originally derived for glasses with quenched 
disorder^ can be extended to quasi-equilibrated super- 
cooled liquids with the result that Lu{(t>) is an inverted 
parabola with a temperature-independent curvature: 

uj{<j)) =uJo ^ . (7) 

Even though combinatorial arguments suggest that 
the parabolic approximation should fail at low 
temperatures,—!^ the high-temperature portion of LL>{if)) 
is supported by simulations of Lennard-Jones (LJ) 
mixturesji^i^ A more stringent test of the Gaussian 
model comes from considering the temperature depen- 
dence of the average basin depth 0y(T) which, in the 
Gaussian model, is predicted to scale linearly with l/T, 
producing a 1/T^ scahng for the configurational heat ca- 
pacity. Most available simulations report deviations from 
this scaling4^i^i^i^ It is currently not fully established 
whether the origin of this deviations should be traced 
back solely to anharmonicity effects^ or to the actual 
failure of the Gaussian approximation although two re- 
cent simulations exploring the low-temperature part of 
the landscape point to the latter explanation42i^ In the 
excitations model, the temperature dependence of the 
average basin energy is qualitatively different for fragile 
and intermediate liquids. In the former case, the basin 
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energy is essentially flat for fragile liquids terminating 
through a discontinuity at the phase transition below 
Tg. For intermediate liquids, 4>v{T) starts to dip as 1/T 
from a high-temperature plateau (as was found for the 
BMLJ liquids"'^'^) and then inflects into an exponential 
temperature dependence recently observed in simulations 
of model network glass-formers,— Note that the excita- 
tions model is the only analytical model we are aware of 
which incorporates both types of temperature scaling in 
one formalism. 

The Gaussian excitations model leads to a non- 
Gaussian P(0) and thus a non-parabolic When 
the non-Gaussian distribution P((/)) is fitted to a Gaus- 
sian function, the result is an approximately linear scal- 
ing of the squared width with temperature CTp y (T) oc 
T. This results makes P{(t>) in Eq. ([6]) a Boltzmann 
distribution, which seems to be more relevant for a 
(quasi)equilibrated supercooled liquid than the non- 
Boltzmann distribution of the Gaussian landscape more 
relevant for systems with quenched disorder. In addi- 
tion, the excitations model gives hyperbolic temperature 
scaling for the heat capacity Cp y oc 1/T. This scaling 
is often observed in the laboratory for c™(T), but here 
we again face the same problem as above concerning the 
connection between excess and configurational heat ca- 
pacity. 

An approximately linear scaling of the width of P{4)) 
with temperature in the excitations model is the result 
of the assumed mean-field, infinite-range attraction be- 
tween the excitations, which is equivalent to assuming a 
Gaussian manifold of real-space excitation energies with 
the variance 2AT. A finite range of interactions between 
the excitations will produce a more complex tempera- 
ture dependence. For instance, a recent exactly solvable 
landscape model of the fluid of dipolar hard spheres^ 
gave a fairly complex temperature scaling of the distri- 
bution variance cFy{T) cx (1 -f (3b)~^ {b is an interac- 
tion parameter). Notice in this regard that a fluid with 
dipolar interactions is an archetypal system in which the 
mean-field approximation is not applicable. The reason 
is two-fold: (i) the average of the potential is zero and 
fluctuations is the first non- vanishing contribution to the 
thermodynamics^ (cf. to LJ forces described reasonably 
by the mean-field van der Waals model) and (ii) the inter- 
action potential is strongly anisotropic. How the model 
should be extended to a finite range of correlations be- 
tween the excitations is not clear now, but the distribu- 
tion width is expected to transform to a temperature- 
independent value for isotropic short-range LJ forces, 
in compliance with the Gaussian model [Eq. ([7])]. Any 
non-Gaussian landscape will generate a temperature- 
dependent width when distribution of inherent energies 
is fitted by a Gaussian. 

In this paper, we use the results of simulations of the 
mSW potential to address some of the challenges listed 
above. We present the results for the landscape thermo- 
dynamics for two values of the tetrahedrality parameter 
A of the mSW potential (see below) and will compare 



the results of the analysis to some other models of glass- 
formers on one hand and to laboratory data for metallic 
glass-formers on the other. Some thermodynamic pa- 
rameters relevant to our discussion are listed in Table [D 
All energies throughout below are in kelvin and entropies 
and heat capacities are in units of k^- Also, we use low- 
case letters for thermodynamic potentials and energies 
per liquid particle, e. g. Spy(T) refers to the configura- 
tional entropy per particle. 

II. LANDSCAPE THERMODYNAMICS 

The properties of the energy landscape for a given in- 
teraction potential can be studied by either a direct cal- 
culation of the enumeration function uj{(j)) or by looking 
at the ensemble averages^i In the first route, w(0) is cal- 
culated by patching together the distribution functions 
at different temperatures once the total and basin free 
energies entering Eq. ^ are knowni^i^i^ The width 
of each individual distribution scales as 1/ Vn with the 
number of particles N, and most simulation data allow a 
Gaussian fit 

P(0) oc e-^(^-^-.WT))V2rp.v(T)^^ 

Here, Tp^viT) is an empirical Gaussian width and the 
stationary point (j)py{T) is the solution of the equation 

/3 (H- df\cj,, T)ld^) = du:{cp)/dcp. (9) 

When f^{(f>,T) = f^{T) is independent of (harmonic 
approximation) and the enumeration function is given by 
the inverted parabola [Eq. ([7])], one gets the hyperbolic 
temperature scaling 

2 

0P,y(r) = 0o-Y (10) 

characteristic of the Gaussian landscape. 

The configurational entropy is determined by the enu- 
meration function taken at the average basin depth 

s'i,y{T)=u;{^Py{T)). (11) 

The configurational entropy is then a part of the thermo- 
dynamic free energy f{T) 

f{T) - 4>py{T) - Ts<i,y{T) + f[4>py{T),T]. (12) 

Once both (j)py{T) and Spy{T) are known, the enumer- 
ation function u;((/)) can be calculated using Eq. pT|) ^ 

III. SIMULATION DETAILS 

We use a model of network liquids which was originally 
introduced by Stillinger and Weber (SW)^'^ for silicon. In 
the SW model, a three-body term is added to the pairwise 
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TABLE I: Thermodynamic parameters (energies in kelvin and entropies in fee units) of mSW and BMLJ models. 



Model 


e" 


Tfn 


rpC b 




To' 


— 010° 


— 0DC^ 




cx h 

■So 


So 


mSW (A = 19) 


25150 


806 


434 


392 


420 


49482 


50300 


1.20 


2.4 


1.0 


BMLJ (p = 1.2) 


120 




35 










0.93 




0.32 



"Lennard- Jones energy of the mSW and BMLJ potentials. 

'From extrapolating Sp{T) to zero for the mSW liquid and from 
the constant-volume dataH for the BMLJ liquid. TJ=- = 487 K 
obtained from s^(r). 

'^From solving the equation s'p'(T) = with excess entropy given 
by Eq. | [26t . 

''From the fit of diffusivity data from MD simulations to the Vogel- 
Fulcher-Tammann equation. 

'^Extrapolated from 4>p{T) obtained from NPT simulations to the 
temperatures at which Sp{T) becomes zero; —(pic = 49265 K 
for the NVT ensemble. 

■'Basin depth of diamond cubic crystal. 

^Top of the enumeration function corresponding to T — > oo limit 
[e.g., ujQ in Eq. Q]. For mSW fluid the number was obtained from 
fitting the simulation data for Sc{T) by a polynomial in 1/T; for 
BMLJ fluid the value listed is from Ref. [l^. 

''Obtained by numerical extrapolation the simulated excess en- 
tropy s'p^(T) of the mSW fluids over the DC crystal to the limit 
T ^ oo. 

'Configurational component of the excitation entropy obtained 
from fitting the conflgurational thermodynamics data from numer- 
ical simulations. 



potential V2{r) to introduce penalty for deviating from 
tetrahedrality 

f(fl2,?'23,?'3l) = V2{ri2) + Au3(ri2,r23,r'3i) (13) 

where 

.2(r)^| MB/ryi)eMr-a)-'lr<a ^^^^ 

with A = 7.049555277, B = 0.602245584, and a = 1.8. 
The three body potential has form 

i'3(?'i2,f 13, 6*213) = exp[7(ri2-a)"^+7(ri3-a)"^]x(3cos6'2 

(15) 

with 7 — 1.2. The potentials are given in reduced units 
(7 = 0.20951 nm and e = 50 kcal/mol (see also Table|I]). 

The original SW model with A = 21 describes sili- 
con, but was modified recently by Molinero et al.^^ by 
decreasing the tetrahedrality parameter A (in contrast to 
its increase attempted earlier by Middleton and Wales^) 
to obtain monoatomic glass-formers. They showed that 
the system crystallizes into diamond cubic (DC) lattice 
for A > 20.25 and into body-centered cubic (BCC) lattice 
for A < 17.5. For intermediate values of A the fluid fails 
to crystallize on the time-scale of computer simulation, 
producing glass-formers. 

Since each fluid, characterized by a given value of A, 
has an equilibrium crystalline phase, this property can be 
used to obtain both the excess and conflgurational ther- 
modynamics for the same system. Two fluids have been 
used in simulations: the original SW model (A = 21) 
and mSW model (A = 19). The results were obtained 
from NVT and NPT MD simulations using the constraint 



method^^ and its modiflcation for the NPT ensemble.— 
Periodic boundary conditions for the cubic cell of 512 par- 
ticles were applied and the time step was about 1.53 fs. 
The temperature has been changed in step-like way with 
50 K (0.002e) per jump. The run length at each given 
temperature varies between 0.76 ns at high temperatures 
to up to 3 ns at the lowest temperature equivalent to 
cooling rates of 65 K/ns and 16 K/ns, respectively. 

IV. RESULTS 

+ 1)V9 

A. Vibrational thermodynamics 

The excitation profiles (j)py{T) and the distribution 
widths Tpy{T) for two SW potentials characterized by 
A = 19 and A = 21 are shown in Fig. U) The average 
basin depth (l)py{T) from simulations was fitted to the 
function 

4>py{T) = a%WXT*W^\{T*f+h%/T*+h%/{T*f 

(16) 

with T* = T/e. The expansion coefficients in Eq. are 
listed in Table [Til The basin depth, shown in Fig. d] rel- 
ative to the equilibrium crystalline (DC) state, changes 
little on this energy scale. The most interesting observa- 
tion is an approximately linear increase of the effective 
width with temperature. This result is inconsistent with 
the Gaussian landscape model [Eq. ([T])]. A part of the 
of the width increase at constant pressure comes from 
thermal expansion (cf. triangles with circles in Fig. d]), 
but there is still an increase of the effective width by a 
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FIG. 4: Average basin energy (a) and the effective Gaussian 
widtli [(b) , Eq. ((8}] vs temperature reduced by the LJ energy 
e (Table [ij. Open points refer to A = 21, closed points in- 
dicate A = 19. Triangles and circles refer to NVT and NPT 
simulations, respectively. The depth of the diamond cubic 
basin ^dc (Table is used for the energy scale. Tm in (a) 
marks the melting temperature and ^ig = (t>p{Tpc) marks the 
depth of the ideal-glass minimum at which the configurational 
entropy Sp{T) becomes zero (A = 19). 



mSW (X=19) 
liquid 



T =0.023 
T* = 0.070 




FIG. 6: Density of states for mSW model (A — 19) in the 
liquid state at constant pressure (a) and DC crystal (b); T* = 
T/e, where e is the LJ energy (Tab.|T]) 
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FIG. 5: Basin free energy for the mSW fluid calculated in the 
harmonic approximation [Eq. (I17|l ] from NPT simulations at 
A = 19, pa^ = 0.51, and varying temperature T* — T/e: 
0.026 (circles), 0.028 (squares), 0.03 (diamonds), 0.036 (left 
triangles), 0.046 (stars). The solid line is a linear regression 
through the points: <j>p,v ~ o-p,v + fcp,v(0/e) with ap — 
-23.63 and bp = -5.80. The NVT data, not shown here, 
are very close to the NPT data with the linear regression 
coefficients ay = —22.11 and by = —5.02. 

factor of nearly 2 even in constant- volume simulations. 
Also shown in Fig. |4] is the energy of the deepest amor- 
phous minimum (j)p{TK), corresponding to the ideal-glass 
transition, measured at the Kauzmann temperature Tk 
at which the configurational entropy Sp (T) becomes zero 
(Table|T|. This minimum lies about 820 K above the crys- 
talline minimum, which can be compared to Stillinger's 
estimate of 460 K for laboratory OTP."'^ The excess val- 
ues are consistent with the requirement that even ideal 
glasses are metastable with respect to the corresponding 
crystals. 

In order to gain insight into the origin of the tempera- 



ture increase of TpyiT) seen in Fig. |4] one needs to sep- 
arate the basin free energy f^{4), T) in Eq. ([6]) from the 
enumeration function. One expects that harmonic ap- 
proximation holds at low temperatures when the basin 
free energy can be obtained by diagonalizing the Hes- 
sian matrix at the local minimum of depth along the 
simulation trajectory 

3Ar-3 

Pf'^y{cl),T) = N-\ ln/3?i^f''')0. (17) 
1=1 

We found, as in previous simulations ,'*^''*^i^° that 
fpy{(j),T) obtained from Eq. (|17p is an approximately 
linear function of (j) (Fig- E]) 

lifp,v{<t>, T) + 31n(r/e) « apy + bpy{c^/e), (18) 

where the linear regression coefficients are listed in the 
caption to Fig. [5l The basins thus become increasingly 
sharp on cooling, both at constant pressure and constant 
volume conditions. The latter observation is an agree- 
ment with the previous constant-volume simulations of 
SPC/E water ^ but in contrast to the opposite trend 
found in constant-volume simulations of BMLJ fluids 

The vibrational density of states (VDOS) at constant 
pressure, used to calculate the harmonic part of the basin 
free energy, is presented on Fig. [51 It was obtained by 
diagonalizing the Hessian of potential energy in inherent 
structures using LAPAC's routine DSYEV,S The VDOS 
of the liquid phase (Fig. [5^) is continuous, with two well- 
defined maxima corresponding to low-frequency longitu- 
dinal and the high-frequency transverse vibrations. The 
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TABLE II: Fitting coefBcients in Eq. (|16|) for the average energy of inherent structures obtained from NPT and NVT MD 
simulations. 

A Ensemble 'Jp.V '^p.V v ^p.V ^p.v- 

A = 19 NPT -2.0013 0.79174 -2.2294 0.0020186 -2.86 x lO"'' 

NVT -1.9744 0.41355 -1.5993 0.001723 -3.03 x lO"'^ 

A = 21 NPT -8.0414 66.40 -264.45 0.25165 -3.82 x 10^^ 

NVT -2.8242 9.5033 -35.642 0.04262 -7.34 x 10"'' 



VDOS shifts to higher frequencies on cooUng, in agree- 
ment with Fig. [51 The crystalhne VDOS was calculated 
for the ideal diamond cubic crystal with = 1728 parti- 
cles. In contrast to the VDOS of the liquid phase, it has 
a discrete spectrum sensitive to the system size. Because 
of this complication, the vibrational thermodynamics of 
the crystal was calculated at different crystal sizes and 
infinite-size results were obtained by linear extrapolation 
of the dependence to iV ^ oo. 

The eigenfrequencies of the Hessian matrix uji depend 
weakly on temperature. This effect is caused by anhar- 
monicity of basin vibrations which tends to soften vibra- 
tional frequencies with increasing temperature. There- 
fore, in order to properly calculate the harmonic free en- 
ergy fpy{T), we used the extrapolation oiuJi{T) to T = 
0. In this case, Eq. p7p with temperature-independent 
frequencies gives the expected value for the harmonic 
part of the basin internal energy 

e'}.y{T) = d{pfi.y)/dp = 3T(1 - l/N). (19) 

Given that the basin free energy f'^{ip,T) is a linear 
function of the basin energy </> (Fig. [5]), the width of the 
Gaussian distribution r(r), obtained by fitting to the 
probability P{(j>) function to Eq. ([5]), is equal to the width 
a obtained by quadratic expansion of the enumeration 
function around the average basin energy 4>p.v{T) (sec- 
ond derivative of /(</>, T) in is zero). This conclusion is 
limited by the neglect of the second derivative of the an- 
harmonic part of the basin free energy /^"''(0, T) which 
we could not extract from our simulations. Since the 
Gaussian landscape precludes temperature dependence 
of the width, the approximately linear increase of the 
width of P{4>) with temperature seen in Fig. 2] can be 
assigned to an increase of the effective width of a non- 
Gaussian enumeration function fitted to a Gaussian. 

The basins of the mSW fluid are anharmonic, as is 
seen from the comparison of the potential energy u{T) 
to the energy 4>py{T) + (3/2)T (Fig. H lower panel). 
The anharmonic potential-energy part, 

u^^^{T) = e(T) - 3T - 4>p,v{T), (20) 

of the internal energy per particle e(T) can be used to 
calculate the anharmonic part of the basin free energy 
according to the thermodynamic equation 

T 

Pfpy{T) = J u^^{T')d(^ljy (21) 
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FIG. 7: Total entropy s (upper panel) and potential energy 
u (lower panel) vs T. The entropy is split into the configura- 
tional entropy s^, harmonic entropy s^, and the anharmonic 
entropy s^^. The potential energy (lower panel, circles) is 
compared to the average basin depth 0y (triangles) and the 
potential energy in the harmonic approximation (jiv + {3/2)T 
(squares) (A = 19, NVT simulation). 



B. Configurational thermodynamics 

Once the harmonic and anharmonic contributions to 
the basin free energy are available from Eqs. (jl7p and 
(|2T1) . the configurational entropy can be calculated by 
subtracting the vibrational (harmonic and anharmonic) 
entropy of basins from the total entropy 

s'pyiT) = s{T) 4y{T) - s|!;{>(T). (22) 

In Eq. (|2ip , the harmonic entropy is calculated from T = 
extrapolated basin frequencies as 

3JV-3 

s%^y{T) = (7V)-i 5] [1 - HPhoj,)] , (23) 
1=1 

and the anharmonic entropy is obtained from Eq. (j2ip . 

Thermodynamic integrationi^i^" was employed to cal- 
culate the total entropy in Eq. ^F^. The excess free 
energy over that of the ideal gas below some reference 
temperature was obtained by integrating the internal 
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energy from simulations: 

A/(T, p) = AfiTr, p) + / e{p, P')dp', (24) 

where Tr was chosen above the critical temperature. The 
value A/(Tr, p) was then obtained by isothermal expan- 
sion to the ideal gas using the equation 



/3,A/(T„p)= / % 
Jo P 



P dp' (l3rP{p') 



1 



(25) 



Finally, the free energy of the ideal gas was added to 
obtain the total free energy of the mSW fluid. 

The excess entropy was calculated from the 
temperature-dependent enthalpies of the liquid and 
the crystal according to the relation 



dT' d 
T' dT' 



ffiiq(r')-i?Dc(r') 



[26) 

where As^p^ = 2.10 (A = 19) is the fusion entropy. The 
enthalpies of the liquid phase, _ffiiq(T), and the DC crys- 
tal, Hdc(T), were fitted from the simulation data to the 
following functions: 



i/iiq(T) = 4>{T) + 3T + C2T^ + C3T\ 
Hdc{T) = -2e + 3T + dsT' + ^^3^^ 



(27) 



where the polynomial coefficients are: C2 = 12.3195e ^, 
C3 = -107.993e-^ da = 0.9352e-\ and = 26.6414e-2. 
The configurational entropy, alternatively to Eq. (P^ . 
can be calculated from the configurational heat capacity 



(T) 



'p.y 



(To)- 



1 d^py{T') 



T' 



dT' 



dT', (28) 



where Tq is some temperature for which Spy (To) is 
known from the thermodynamic integration to the ideal 
gas. The two thermodynamic routes give identical re- 
sults. 

The splitting of the total liquid entropy into vibra- 
tional and configurational parts is shown in the upper 
panel of Fig. [71 In addition, the configurational en- 
tropy is compared in Fig. [H] to the excess entropy over 
the thermodynamically stable DC crystal in the temper- 
ature range between the melting temperature and the 
lowest temperature accessible to simulations. Also shown 
are the harmonic and anharmonic parts of the excess vi- 
brational entropy. The vibrational entropy makes about 
half of the overall excess entropy, in a general accord 
with experimental evidence obtained for laboratory glass- 
formers<22i22iiHiS The situation is somewhat similar with 
the excess heat capacity close to the melting point where 
the anharmonic vibrational part is responsible for ap- 
proximately half of the excess heat capacity (harmonic 
excess heat capacity is identically zero). However, the 



2.5 
2 

1.5 
1 

0.5 


2.5 

2 

1.5 









ex / 

Sp (a) 


- IGfits 


anh 


vib 















\ C 






T 

ex m - 


IG fit - - 







550 600 650 700 750 800 
T/K 

FIG. 8; (a) Excess entropy of the mSW fluid over its DC crys- 
tal (s*^) and its configurational (spy) and vibrational (s^''') 
components. The vibrational excess entropy is split into har- 
monic (sp,v) 8.nd anharmonic (sp^^) parts. Configurational 
entropy at constant volume is shown by the dash-dotted line, 
(b) The splitting of the excess heat capacity into configura- 
tional (cp) and anharmonic vibrational (cp"*^) contributions. 
The dashed lines in (a) and (b) refer to fits of the configura- 
tional thermodynamics to the excitations (IG) model.— The 
fit parameters {eo,A,so} are: {1675 K, 787 K, 1.0}. NVT 
simulations were done at density pa^ = 0.51 and the density 
changes from 0.524 at Tm to 0.528 in NPT simulations. 



fraction of the anharmonic heat capacity drops down to 
about 15% with lowering temperature (Fig.[5jD). 

We have used the data for the configurational thermo- 
dynamics from simulations to fit them to the Gaussian 
excitations (IG) model?^ The model does not anticipate 
anharmonicity playing a major role in the excess ther- 
modynamics below the melting temperature. Our focus 
is therefore limited to configurational thermodynamics 
only. As is shown in Fig. [U the IG models can be suc- 
cessfully fitted to the temperature dependence of the con- 
figurational entropy and to the high-temperature portion 
of the heat capacity. The model, however fails to repro- 
duce the sharp rise of the heat capacity at the lowest 
temperatures accessible to simulations. 



C. Dynamics 

The diffusivity data from simulations (A 19) are 
shown by points in Fig. [9l These results are fitted to 
the AG relation [Eq. ([T|)] and to the dynamic equation of 
the Gaussian excitations model [Eq. The configura- 
tional heat capacity from our simulations is used in the 
fit, in contrast to the previous application of Eq. ([5]) (Ref . 
[33 ) where experimental dielectric relaxation daia^ were 
fitted to Eq. ([5]) using the excess heat capacity from the 
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FIG. 9: Fit of the diffusivity of the mSW (A = 19) fluid 
from simulations (hatched squares) to the AG theory [dashed 
line, Eq. ^] and to the excitations model [solid line, Eq. ((S))]. 
The dashed-dotted line refers to the Vogel-Fulcher-Tammann 
(VFT) equation with the VFT temperature To equal to 303 
K (Table H)). The dashed and sohd lines are indistinguishable 
on the scale of the plot. 



laboratory experiment. However, for the mSW model, 
c™{T) and Cp{T) are ofF-set by almost a constant shift 
of anharmonic heat capacity, and the use of either of 
the two to fit diffusivity gives comparable results. The 
dashed line, almost indistinguishable from the solid line 
in Fig. [31 indicates the AG relation. The Vogel-Fulcher- 
Tammann (VFT) equation (dash-dotted line in Fig. [3]) 
gives a less satisfactory fit. In terms of fitting the re- 
laxation data, the AG relation is superior to both the 
Gaussian excitations model and the VFT equation since 
it involves one fitting parameter less. 



V. DISCUSSION 

In discussing these results we must first recognize the 
sort of frustrations that are likely to accompany any ef- 
fort to resolve the key problems of the glass transition 
by studying non-crystallizing systems using MD meth- 
ods. Despite the four orders of magnitude in diffusivity 
that we have studied (Fig. [9]) we have barely reached the 
onset of the "low temperature domain" (from the Stickel 
temperature Tb down to Tg) in which the Adam-Gibbs 
equation has been tested experimentally using the excess 
entropy. It is in this domain that experiments show lin- 
ear relations between logD and {Ts°p)~^ predicted by 
the AG equation. Thus when, in Fig. [TOl we plot our 
\ogD values against the alternative quantities {Tsp)^^ 
and (Ts^)~^, and observe that the first is linear and 
the second is not, we are not able to relate the break in 
the (Ts^)^^ plot to the breakdown of the AG correla- 
tion at Tb in the analysis of Richert and Angell,^'^ nor 
to the other crossovers (Stokes-Einstein equation break- 
down etc.) observed in the experimental plots at Tc j^^i^° 
The break in our logZ? vs (Ts^)~^ plot occurs at a 
quite different (much higher) temperatures, where D is 
only 10~^ cm^s~^, and thus must have a different ori- 
gin. Whether or not there is a further break in the plot 
of \ogD vs (Tsf)-^ (or of logD vs (Ts'p)-^ for that 
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FIG. 10: (a) AG plot for diffusivity using excess (triangles) 
and configurational (circles) entropy of the mSW fluid with 
A = 19. (b) s5?(r) vs s'piT) for the same system. The dashed 
lines represent the values of corresponding parameters at the 
melting temperature. 



matter) occurring at Tj, cannot be told from the present 
work, nor from any previous study. 

This difference between the viscosity domains explored 
in MD simulations of glass-formers and the low temper- 
ature domain near Tg, where so many laboratory stud- 
ies are carried out, is not given adequate attention in 
most of the discussions of "glassy dynamics" simulation 
results. Notwithstanding the success of well know phe- 
nomenological models in describing the glass transition 
as is observed in simulation, the fact that one is work- 
ing above the much discussed crossover in the MD case 
and below it in the experimental range of AG equation 
testing, cannot be escaped. The best that can be done is 
to compare the slopes of the two plots with those found in 
the most relevant experiments, but little can be gained 
thereby without a better theory for the AG equation. 
Here wc will further discuss the break in our Fig. [10] plot 
for diffusivities and the non-Arrhenius character of the 
diffusivities, and will then seek to reconcile what seems 
to be a conflict in dynamic and thermodynamic signa- 
tures of fragility in the present system. This will provide 
us with the opportunity to make a (rare) comparison of 
thermodynamic behavior for different potential models 
of simulated glass-formers. 

The most complete studies of glass-former diffusivity 
available are for the cases of OTP;^ Si02r^ and some 
of the bulk metallic glasses (BMG)^iH in which data 
cover the range of diffusivity from water-like values down 
to those characteristic of liquids at their glass transition 
temperatures (10~^^ cm^s~^). Only in the case of OTP 
is the variation of the excess entropy in the same tem- 
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FIG. 11: Diffusivity of model and laboratory glass-formers 
vs Tr/T, where Tr is the temperature at which diffusivity 
is equal to 10~^ cm^s~^ (the end of the Arrhenius region). 
Model fluids: larger component of the BMLJ liquidi^ and 
mSW (A = 19) fluid. Laboratory Uquids: ^'^Co (triangles)^ 
and Ni (pluses)^ tracers in Pd43Cu27NiioP2o melt and OTP 
(squares).— Stars refer to the diffusivity of ^''Co calculated 
from the Stokes-Einstein equationi^ 



perature range, relative to that of the crystal, properly 
known. Excess entropies relative to a mixture of crystals, 
are known for some of the bulk metallic glass-formers. 

The division of the excess entropy (and heat capacity) 
of the supercooled liquid into vibrational and configura- 
tional components was suggested in Goldstein's original 
analysis, where it was found that in the case of OTP 
almost 50% of the excess entropy was vibrational in char- 
acter. Goldstein's finding has recently been confirmed by 
measurements of Wang and Richerti^ However OTP is 
fragile in character and the behavior of the VDOS, in 
Lewis- Wahnstrom model, is unlike that of the present 
system so comparison of our findings with those for OTP 
may not be appropriate. The bulk metallic glasses, by 
contrast, are more similar to the present system in their 
VDOS behavior (from neutron scattering studies of their 
quenched and annealed states,^ but remember the ob- 
servations were all made at fictive temperatures near Tg) 
and prove to be relatively strong in their kinetics Like 
the present system, their diffusivities exhibit a strong Ar- 
rhenius plot curvature in the temperature range acces- 
sible to computer simulation (thus appearing fragile in 
this range) in much the same way as do classical network 
glasses, BeF2 and Si02, at high temperatures. Thus the 
behavior of our mSW system might be better compared 
with that of the BMG systems studied by Chathoth et 
ali^ and reviewed by Faupel et ali^ 

In Fig. [TI] we compare the diffusivities from Fig. [H] 
with those of various components of the bulk glass-former 
Pd43NiioCui3P2o from Refs. [sl and HI after scahng by 
the temperature at which each system exhibits D = lO"'"^ 
cm^s~^ (near where the deviation from Arrhenius behav- 
ior first becomes obvious). On the larger temperature 
scale the Ni diffusivity, like the ^^Co diffusivity of Ref. 
l63l and the viscosity of Ref. [63|, all return to Arrhenius 
behavior with a larger slope, and the behavior appears 
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FIG. 12: Excitation profiles {(t>viT)) of the mSW liquid (dia- 
monds, A = 19) vs temperature measured on the Kauzmann 
temperature scale; T^{V = Const) is from extrapolating 
Sy{T) to zero. Solid line is the excitation profile obtained 
from the flt of IG model to the constant- volume configura- 
tional entropy of the mSW liquid obtained by thermodynamic 
integration from the ideal gas reference state, and horizontal 
dashed line marks the position of the ideal glass basin depth. 
Open circles are data from the slowest effective cooling rate 
set of inherent structure energies from the study of Sastry et 
al on BMLJ at pa'' = 1.213~ scaled using the Kauzmann 
temperature = 34.35.— The level at low T in this case 
is not the ground state but the glassy state frozen in at this 
cooling rate. The sharper profile for the mSW model implies 
a more fragile liquid than the BMLJ model, as is also seen in 
the diffusivity data in Fig. 1111 



non-fragile, approximately like glycerol. Thus the strong 
curvature which lead Molinero et ali^ to conclude that 
mSW is very fragile, does not necessarily continue. This 
would rationalize what otherwise is a problem raised by 
the excitation profile of Fig. 0]- which is not of the form 
expected for a very fragile liquid according to the IG 
model. This point is illustrated in Fig. [T^ where the ex- 
citation profiles of the BMLJ and mSW liquids are scaled 
onto the same plot by the use of their Kauzmann tem- 
peratures (r|. = 34.35 Ki3 for BMLJ and Tk = 434 
K for mSW, Table |T|. Unlike the more fragile casesjiii^^ 
which develop an S-shaped profile and exhibit phase tran- 
sitions (not unlike that of silicon itself), these profiles al- 
ways have positive slopes. Consistent with the difference 
in their diffusivity behavior seen in Fig. Ill[ the profile 
for mSW is sharper than that of the less fragile BMLJ. 
Figures [TT] and [12] together provide the best evidence to 
date of the surprisingly non-fragile behavior of the much 
studied BMLJ system. 

A way of inducing fragile character in an atomic sys- 
tem spherically symmetric potential is, according to Sas- 
try, to increase the density of BMLJ. This was shown to 
increase the slope of the AG plot, Fig. [TUb . in the same 
way that is seen when the experimental data for OTP 
are added to the plot. We demonstrate this in Fig. [T31 
To show consistency, we include, in Fig. [T31 the data for 
bulk metallic glasses, using the excess entropy data of 
the BMG reported by Kuno et ali^ Although these data 
refer to the melting of ternary eutectic composition, we 
consider that the fusion enthalpy used in the entropy as- 
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FIG. 13: Comparison of AG plots of laboratory glass-formers 
to the simulation results of mSW (A = 19). Laboratory glass- 
formers: OTP (crosses^); Ni (diamonds^) and ^^Co (closed 
squares--) tracers in Pd43Cu27NiioP20 melt. Open squares 
refer to the diffusivity of the same melt calculated from vis- 
cosity by using the Stokes-Einstein equation.— 



sessment to be valid (i.e. to contain negligible non-ideal 
mixing enthalpy), because it has been shown that the 
crystals that fuse at the eutectic are already binary com- 
pounds. The slope of the plot for the BMG diffusivities 
for Ni and Co, which are decoupled from the viscosity, 
is less than for mSW using excess entropy, but the slope 
for the viscosity-based data of Fig. [TT] is essentially the 
same. The variations in slope in Fig. 13, however, are 
not well accounted for. Here we recall that, in Ref. [3^ . 
such plots could be reduced to an all common slope using 
Eq. (O of this paper. 

The analysis of the thermodynamic and relaxation 
data of laboratory glass-formers has allowed us to clas- 
sify those liquids according to thermodynamic parame- 
ters of their configurational excitationsiiii^ The excita- 
tions model describes the thermodynamics of a super- 
cooled liquid as an ideal gas of excitations each carrying 
excitation energy e and entropy sq. The energy of exci- 
tations belongs to a Gaussian manifold with the average 
excitation energy eg ~ 2x(T)A {x(T) is the population of 
the excited state) and the variance 2TA. This represen- 
tation is in fact equivalent to an ensemble of excitations 
with mean-field attractions. When the laboratory data 
for fragile liquids are fitted to the IG model, the excita- 
tion parameters show universality when the energies are 
scaled with the Kauzmann temperature. In these reduced 
units, the excitation entropy sq becomes the only rele- 
vant parameter determining fragility as shown in Fig. [T3] 
where the data are plotted against the steepness fragility 
index;^ 



A 



2Tk Tk 



So- 



(29) 



This universality does not hold for liquids of interme- 
diate fragility (intermediate liquids). Noteworthy is a 
much smaller width parameter A for intermediate liquids 
compared to fragile liquids (Fig. [14)) . 

The current simulations of the mSW liquid generally 
support the basic assumptions incorporated in the de- 
velopment of the excitations model of low-temperature 
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FIG. 14: Excitation energy eo, the trapping energy param- 
eter A, and the excitation entropy so of some laboratory 
glass-formers. The energy parameters are scaled with the 
experimental Kauzmann temperature Tk- The parameters 
obtained by fitting the excitation IG model to experimen- 
tal excess thermodynamic data are plotted against the steep- 
ness fragility parameter m.— The fragility number m ~ 60 
separates fragile from intermediate liquids (dash-dotted line). 
The numbers in the plot indicate: toluene (1), D,L-propene 
carbonate (2), OTP (3), 2-methyltetrahydrofurane (4), salol 
(5), 3-bromopentane (6), glycerol (7), selenium (8), and n- 
propanol (9). For selenim, the steepness index of m = 40 
calculated in Ref. (63 was used. 



glass-formerSfiii^^ although, as for the comparison to lab- 
oratory data, the MD simulations probe the temperature 
range much higher than the one anticipated in the theory 
development {Tg < T <Tb). Nevertheless, some conclu- 
sions can be drawn. The IG model neglects the temper- 
ature dependence of the excitation entropy and thus the 
anharmonicity effects. Anharmonic vibrations are signifi- 
cant in the potential energy landscape of the mSW model 
at high temperatures, but their contribution to the excess 
heat capacity drops to about ~ 15 % in the lowest por- 
tion of the temperature range studied here. Therefore, 
the assumption of the temperature-independent entropy 
of elementary excitations might be a good first-order ap- 
proximation in applications of the model to data at low 
temperatures close to Tg. 

The most significant ingredient of the excitations 
model requiring testing on available data is the antic- 
ipated temperature dependence of the effective Gaus- 
sian width of the distribution of basin energies. The 
temperature-dependent width is introduced into the 
model to account for an effectively non-Gaussian land- 
scape probed by the system when exploring deeper basins 
with lowering temperature. An exact non-Gaussian land- 
scape model recently developed by one of us^ has taught 
us that this temperature dependence can be quite com- 
plex such that a linear scaling probably applies only to a 
limited range of temperatures. Even though, in the range 
of temperatures studied by our simulations, the effective 
width scales linearly with temperature for both constant- 
volume and constant-pressure ensembles (Fig. ID). Unfor- 
tunately, the present simulations do not allow access to 
the lower portion of the landscape, and deviations from 



gaussianity observed in some recent simulation; 



,49.51 



are 
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not pronounced in the enumeration function a;(^). 0082535 (C. A. A.). We are grateful to Srikanth Sastry 

for sharing his experience with numerical calculations. 

Acknowledgments 

This work was supported by the NSF through the 
grants CHE-0616646 (V. K. and D. V. M.) and DMR- 



1 K. L. Ngai, J. Non-Cryst. Sol. 275, 7 (2000). 

^ C. A. Angell, Science 267, 1924 (1995). 

^ M. Goldstein, J. Chem. Phys. 51, 3728 (1969). 

* G. Adam and J. H. Gibbs, J. Chcm. Phys. 43, 139 (1965). 
^ X. Xia and P. G. Wolyncs, Proc. Nat. Acad. Sci. 97, 2990 

(2000) . 

X. Xia and P. G. Wolynes, J. Phys. Chem. B 105, 6570 

(2001) . 

^ V. Lubchenko and P. G. Wolynes, J. Chem. Phys. 121, 

2852 (2004). 
® H. Basslcr, Phys. Rev. Lett. 58, 767 (1987). 
V. I. Arhipov and H. Basslcr, J. Phys. Chem. 98, 662 
(1994). 

1° J. C. Dyre, Phys. Rev. B 51, 12276 (1995). 

" D. V. Matyushov and C. A. Angell, J. Chem. Phys. 123, 

034506 (2005). 

C. Dallc-Fcrricr, C. Thibierge, C. Abla-Simionesco, 
L. Bcrthicr, G. Biroli, J.-P. Bouchaud, F. Ladiou, 

D. L'Hote, and G. Tarjus, Phys. Rev. E 76, 041510 (2007). 
S. Sastry Nature 409, 164 (2001). 

" I. Saika-Voivod, P. H. Poole, and F. Sciortino, Nature 412, 
514 (2001). 

S. Mossa, E. LaNavc, H. E. Stanley, C. Donati, 
F. Sciortino, and P. Tartagha, Phys. Rev. E 65, 041205 

(2002) . 

N. Giovambattista, S. V. Buldyrev, F. W. Starr, and H. E. 

Stanley, Phys. Rev. Lett. 90, 085506 (2003). 

I. Saika-Voivod, F. Sciortino, and P. H. Poole, Phys. Rev. 

E 69, 041503 (2004). 
^® Y. Gebremichael, M. Vogcl, M. N. J. Borgroth, F. W. Starr, 

and S. C. Glotzer, J. Phys. Chem. B 109, 15068 (2005). 

C. A. Angell, J. Non-Cryst. Solids 131-133, 13 (1991). 
2° R. Richert and A. C. Angell, J. Chem. Phys. 108, 9016 

(1998). 

21 F. H. Stillinger and T. A. Weber, Phys. Rev. A 25, 978 

(1982). 

22 M. J. Goldstein, J. Chem. Phys. 64, 4767 (1976). 

2^ W. A. Phillips, U. Buchenau, N. Niicker, A.-J. Dianoux, 

and W. Retry Phys. Rev. Lett. 63, 2381 (1989). 
2" S. Corczzi, L. Gomez, and D. Fioretto, Eur. Phys. J. E 14, 

143 (2004). 

25 L.-M. Martinez and C. A. Angell, Nature 410, 663 (2001). 
2" C. A. Angell and S. Borick, J. Non-Crystal. Solids 307- 

310, 393 (2002). 
2^ F. W. Starr, C. A. Angell, and H. E. Stanley Physica A 

323, 51 (2003). 
2* G. P. .Johari, ,J. Chcm. Phys. 126, 114901 (2007). 
29 L-R. Lu, G. P. Gorier, and R. Willncckcr, Appl. Phys. Lett. 

80, 4534 (2002). 
^° C. T. Moynihan and C. A. Angell, J. Non-Crystal. Sol. 

274, 131 (2000). 
^1 W. Kob and H. C. Andersen, Phys. Rev. Lett. 73, 1376 



(1994). 

^2 D. V. Matyushov and C. A. Angell, J. Chem. Phys. 126, 
094501 (2007). 

3^ C. A. Angell, Y. Yuc, L.-M. Wang, ,J. R. D. Copley 
S. Borick, and S. Mossa, ,J. Phys.: Condons. Matter 15, 
S1051 (2003). 

L.-M. Wang and R. Richert, Phys. Rev. Lett. 99, 185701 
(2007). 

^5 L. ,J. Lewis and G. Wahnstrom, Phys. Rev. E 50, 3865 

(1994). 

V. Molinoro, S. Sastry, and C. A. Angell, Phys. Rev. Lett. 
97, 075701 (2006). 
^'^ A. Scala, F. W. Starr, E. LaNave, F. Sciortino, and H. E. 
Stanley Nature 406, 166 (2000). 

P. K. Ghorai and D. V. Matyushov, J. Phys. Chem. B 110, 

1866 (2006). 

39 L. Xu, P. Kumar, S. V. Buldyrev, S.-H. Chen, P. H. Poole, 
F. Sciortino, and H. E. Stanley, Proc. Natl. Acad. Sci. 102, 

16558 (2005). 

^° C. A. Angeh and J. Wong, J. Chem. Phys. 53, 2053 (1970). 

*i F. Stillinger, J. Phys. Chem. B 102, 2807 (1998). 

'*2 F. W. Starr, S. Sastry E. LaNave, A. Scala, H. Eu- 
gene Stanley, and F. Sciortino, Phys. Rev. E 63, 041201 
(2001). 

"3 B. Derrida, Phys. Rev. B 24, 2613 (1981). 

K. H. Fischer and J. A. Hertz, Spin Glasses (Cambridge 
University Press, 1999). 

F. H. Stillinger, J. Chem. Phys. 88, 7818 (1988). 

M. S. Shell and P. G. Debenedetti, Phys. Rev. E 69, 051102 

(2004). 

S. Biichner and A. Hcucr, Phys. Rev. E 60, 6507 (1999). 

F. Sciortino, J. Stat. Mechanics p. 05015 (2005). 
*9 A. J. Moreno, I. Saika-Voivod, E. Zaccarelli, E. LaNave, 

S. V. Buldyrev, P. Tartaglia, and F. Sciortino, J. Chem. 

Phys. 124, 204509 (2006). 
5° S. Sastry J. Phys.: Condens. Matter 12, 6515 (2000). 
" D. V. Matyushov, Phys. Rev. E 76, 011511 (2007). 

52 M. A. Osipov, P. I. C. Teixeira, and M. M. T. da Gama, 
J. Phys. A: Math. Gen. 30, 1953 (1997). 

53 F. H. Stillinger and T. A. Weber, Phys. Rev. B 31, 5262 
(1985). 

T. F. Middleton and D. J. Wales, Phys. Rev. B 64, 024205 
(2001). 

55 M. P. Allen and D. J. Tildesley, Computer Simulation of 
Liquids (Clarendon Press, Oxford, 1996). 

5« D. Brown and J. H. R. Clarke, Mol. Phys. 51, 1243 (1984). 

5^" E. Anderson, Z. Bai, C. Bischof, S. Blackford, ,J. Demmel, 
J. Dongarra, J. D. Croz, A. Grcciibauni, S. Harnmarling, 
A. McKcnncy, ct al., LAPACK users' guide (Philadelphia 
: Society for Industrial and Applied Mathematics, 1992). 

5* C. A. Angell, J. Phys.: Condens. Matter 16, S5153 (2004). 

59 M. Mapes, S. Swallen, and M. Ediger, J. Phys. Chem. B 



14 



110, 507 (2006). 

E. Rossler, Phys. Rev. Lett. 65, 1595 (1990). 

N. Giovambattista, C. A. Angell, F. Sciortino, and H. E. 
Stanley, Phys. Rev. Lett. 93, 047801 (2004). 
G. Brebcc, R. Seguin, C. Sella, J. Bevenot, and J. C. Mar- 
tin, Acta Metallurgica 28, 327 (1980). 
A. Bartsch, K. Ratzkc, F. Faupel, and A. Meyer, Appl. 
Phys. Lett. 89, 121917 (2006). 

F. Faupcl, W. Frank, M.-P. Macht, H. Mchror, V. Naun- 
dorf, K. Ratzkc, H. R. Schobcr, S. K. Sharma, and H. Tc- 
ichlcr. Rev. Mod. Phys. 75, 237 (2003). 

A. Meyer, J. Wuttke, W. Petry, A. Peker, R. Bormann, 



G. Coddens, L. Kranich, O. G. Randl, and H. Schober, 
Phys. Rev. B 53, 12107 (1996). 

S. M. Chathoth, A. Meyer, M. M. Koza, and F. Juranyi, 
Appl. Phys. Lett. 85, 4881 (2004). 

S. Sastry, P. G. Debenedetti, and F. H. Stllinger, Nature 
393, 554 (1998). 

M. Kuno, L. A. Shadowspeaker, J. Schroers, and R. Busch, 
Mat.Res. Soc. Proc. 806, 227 (2004). 
L.-M. Wang, C. A. Angell, and R. Richert, J. Chem. Phys. 
125, 074505 (2006). 



